! HYDRO 2D PSEUDOSPECTRAL CODE

! Solves vorticity equation (Navier-Stokes) with optional forcing (Type I)
! \partial_t w = f - \div (w {\bf v}) + \nu \lapl w

! USES:
! COMPLEX ARRAYS, FFTW

! 2nd ORDER RK TIME INTEGRATION
! Full 2/3 De-aliasing or Pseudospectral scheme chosen by external flag

! This version based on previous MHD2D code by Pablo MAR 2001
! Pablo Dmitruk SEP 2001
! <pablo@bartol.udel.edu>

!=====================================================================

  MODULE size_params
    INTEGER ::   nx, ny
    INTEGER ::   nx2, nx1
    INTEGER ::   n2d,n2dh
    SAVE
  END MODULE size_params
!=====================================================================

  MODULE hd_params
      REAL :: xnu
      REAL :: f0, kf0, kf01, kf02
      REAL :: units
      SAVE
  END MODULE hd_params
!=====================================================================

  MODULE fftw_plans
!   For FFTw to work correctly, the size of these integers must
!   match the size of pointers in the language c.
      INTEGER, PARAMETER   ::   ikind = 8       !.. For alpha/64-bit machines
!     INTEGER, PARAMETER   ::   ikind = 4       !.. For most 32 bit machines
!     INTEGER, PARAMETER   ::   ikind = 7       !.. For g77 use (eg, PSR)
      INTEGER( kind=ikind) ::   plan1Df, plan1Db, plan2DRC, plan2DCR
!     INTEGER              ::   plan1Df, plan1Db, plan2DRC, plan2DCR
      SAVE
  END MODULE fftw_plans
!=====================================================================

	PROGRAM hd2d

        USE size_params
        USE hd_params
        USE fftw_plans
        IMPLICIT NONE
        INCLUDE 'fftw_include.h'
        
        COMPLEX, ALLOCATABLE, DIMENSION (:,:)    ::  wk,wh
        REAL, ALLOCATABLE, DIMENSION (:,:)       ::  k2
        REAL, ALLOCATABLE, DIMENSION (:)         ::  kx,ky
        REAL, ALLOCATABLE, DIMENSION (:)         ::  eku

        REAL    :: k2lim
        INTEGER :: istart,nstep,nstg,nstdat,nstk,icnt
        INTEGER :: idale, itime, ii, irk, iwk
        REAL    :: dt, dth, time0, time, ynorm, tiny
        CHARACTER*4  buffe4
        CHARACTER*3  buffe3
        CHARACTER*2  buffe2
        CHARACTER*1  buffe1
        CHARACTER*4  strit 


        REAL    :: fct
        INTEGER :: i,j,ny2

        REAL    :: u0

!       Read input parameters
        OPEN(15,file='hd.inp')
          READ(15,*) istart
          READ(15,*) nx
          READ(15,*) ny
          READ(15,*) xnu
          READ(15,*) dt   
          READ(15,*) nstep
          READ(15,*) nstg
          READ(15,*) nstdat
          READ(15,*) nstk
          READ(15,*) f0
          READ(15,*) kf0
          READ(15,*) units
          READ(15,*) idale
          READ(15,*) irk
          READ(15,*) iwk
        CLOSE(15)
  
        dth=dt/2.
        nx2=nx/2
        nx1=nx2+1

        ALLOCATE( wk(nx1,ny), wh(nx1,ny) )
        ALLOCATE( k2(nx1,ny), kx(nx1), ky(ny) )
        ALLOCATE( eku(nx2) )

        CALL set_plan2D (nx,ny)                ! FFTW plans
        ynorm = 1.0 / (float(nx)*float(ny))    ! FFT normalization factor 
        n2d = ny * (nx + 2)
        n2dh = n2d / 2
        units = ynorm * units
        tiny = 1.e-6
        kf01 = (kf0-1)*(kf0-1)
        kf02 = kf0*kf0

!------------------------------------------------------------
!       Set truncation circle in k2lim:
!       Use this value for removal of double alias errors only:
!       k2lim = (2.0/9.0) * FLOAT( nx*nx)
!       Intermediate Cutoff:
!       k2lim = (1.0/8.0) * FLOAT( nx*nx)
!       Eliminates all aliasing errors:
!       k2lim = (1.0/16.0) * FLOAT( nx*nx)
!       Aliasing removal limit from RR Pablo's code
        k2lim = (1.0/9.0) * FLOAT( nx*nx)
!------------------------------------------------------------

        CALL init_kays (kx, ky, k2)



        IF (istart.eq.0) THEN
			time0 = 0.
			wk = 0.
			u0 = 1.
			!Generates initial stream function 'psi' in real space
			CALL initial_state_real(wk,u0)
			!Transform initial 'psi' to fourier space 
			CALL rfftwnd_f77_real_to_complex(plan2DRC,1,wk,1,n2d,0,0,0) 
			!Obtains initial vorticity in fourier space
			wk = ynorm*k2*wk

			!Write initial state
			IF (iwk.eq.0) THEN
				wh=wk
				CALL rfftwnd_f77_complex_to_real(plan2DCR,1,wh,1,n2dh,0,0,0)            
				OPEN(1,file='w_0000',form='unformatted')
				CALL write_xarr (wh,time0)
				CLOSE(1)
			ELSE
				OPEN(1,file='wk_0000',form='unformatted')
				WRITE(1) time0,wk
				CLOSE(1) 
			ENDIF

		ELSE
			!Restart from previous state
			IF (irk.eq.0) THEN
				OPEN(1,file='w_0000',form='unformatted')
				CALL read_xarr(wh,time0)
				CLOSE(1)
				CALL rfftwnd_f77_real_to_complex(plan2DRC,1,wh,1,n2d,0,0,0) 
				wk=ynorm*wh
			ELSE
				OPEN(1,file='wk_0000',form='unformatted')
				READ(1) time0,wk
				CLOSE(1)
			ENDIF
        ENDIF

		!ale
		wh=wk
		CALL rfftwnd_f77_complex_to_real(plan2DCR,1,wh,1,n2dh,0,0,0)            
		OPEN(1,file='w_0000')
		CALL write_xarr (wh,time0)
		CLOSE(1)
		!ale

!       Impose reality conditions
        CALL preal2D (wk)

!       Write initial global values (if beginning from t=0)
        IF (time0.lt.tiny) THEN
         OPEN(2,file='globals.out')
         WRITE(2,*) nx,ny,nstep,nstg,dt,xnu,f0,kf0
         CLOSE(2)
         CALL globals(wk,k2,time0)
!       Write initial spectra
         CALL spec2D (wk,k2,eku)
!          eku = kinetic  energy spectrum (contains 1/2 factor)
           OPEN(1, file='spec_0000.out')
           WRITE(1,*) time0,nx2
			do ii = 1, nx2 !ale
				write(1,*) ii, nx2
			end do
           !WRITE(1,*) (eku(ii), ii=1,nx2) !ale
           CLOSE(1)
        ENDIF

!       Forcing normalization 
        ny2 = ny/2
        fct=0.
        do j=1,ny
         do i=1,nx1
!        Count only the energy of the mode if ( kx > 0 .or. ky > 0 )
         if ( i.gt.1 .or. (j.gt.1 .and. j.le.ny2) ) then
         if ((k2(i,j).gt.kf01).and.(k2(i,j).lt.kf02)) &
            fct = fct + f0*f0 + f0*f0
         endif
         end do
        end do
        fct = 2.*fct
!       print*,'forcing = ',fct
        
!=====================================================================
!                                                                   
!                  TIME STEPPING STAGE BEGINS                       
!                                                                   
!=====================================================================
!                 2nd ORDER RUNGE-KUTTA SCHEME                      
!=====================================================================

      time=time0
      DO icnt=1,nstep
!
!     HALF STEP
!     =========                                                        
      CALL evaluate(wk,wh,kx,ky,k2,dth,1)
!              wh=wk+dth*dw(wk)

         CALL preal2D (wh)

!        Isotropic truncation:
         IF (idale .eq. 1) THEN
            CALL circ2D (wh,k2,k2lim)
         ENDIF

!     FULL STEP
!     =========                                                        
      CALL evaluate(wh,wk,kx,ky,k2,dt,0)
!              wk=wk+dt*dw(wh)

         CALL preal2D (wk)

!        Isotropic truncation:
         IF (idale .eq. 1) THEN
            CALL circ2D (wk,k2,k2lim)
         ENDIF

      time=time0+icnt*dt

!     Write globals
      IF (mod(icnt,nstg) .eq. 0) CALL globals(wk,k2,time)

!     Write images
      IF (mod(icnt,nstdat) .eq. 0) THEN
        itime=int((time+1.e-4)/(nstdat*dt))
        if (itime.lt.10) then
           write(unit=buffe1,fmt='(i1)') itime
           read(unit=buffe1,fmt='(a1)') buffe1
           strit='000'//buffe1
        else
           if (itime.lt.100) then
              write(unit=buffe2,fmt='(i2)') itime
              read(unit=buffe2,fmt='(a2)') buffe2
              strit='00'//buffe2
           else
              if (itime.lt.1000) then
                 write(unit=buffe3,fmt='(i3)') itime
                 read(unit=buffe3,fmt='(a3)') buffe3
                 strit='0'//buffe3
              else
                 write(unit=buffe4,fmt='(i4)') itime
                 read(unit=buffe4,fmt='(a4)') strit
              endif
           endif
        endif
        IF (iwk.eq.0) THEN
			wh=wk
			CALL rfftwnd_f77_complex_to_real(plan2DCR,1,wh,1,n2dh,0,0,0)            
	!         OPEN(1,file='w_'//strit,form='unformatted') !ale
			OPEN(1,file='w_'//strit) !ale
			CALL write_xarr (wh,time)
			CLOSE(1)
        ELSE
			OPEN(1,file='wk_'//strit,form='unformatted')
			WRITE(1) time,wk
			CLOSE(1)
        ENDIF
      ENDIF

!     Write spectra
      IF (mod(icnt,nstk).eq.0) THEN
           CALL spec2D (wk,k2,eku)
!          eku = kinetic  energy spectrum (contains 1/2 factor)
           itime=int((time+1.e-4)/(nstk*dt))
        if (itime.lt.10) then
           write(unit=buffe1,fmt='(i1)') itime
           read(unit=buffe1,fmt='(a1)') buffe1
           strit='000'//buffe1
        else
           if (itime.lt.100) then
              write(unit=buffe2,fmt='(i2)') itime
              read(unit=buffe2,fmt='(a2)') buffe2
              strit='00'//buffe2
           else
              if (itime.lt.1000) then
                 write(unit=buffe3,fmt='(i3)') itime
                 read(unit=buffe3,fmt='(a3)') buffe3
                 strit='0'//buffe3
              else
                 write(unit=buffe4,fmt='(i4)') itime
                 read(unit=buffe4,fmt='(a4)') strit
              endif
           endif
        endif
           OPEN(1, file='spec_'//strit//'.out')
           WRITE(1,*) time,nx2
			do ii = 1, nx2 !ale
				write(1,*) ii, nx2
			end do
           !WRITE(1,*) (eku(ii), ii=1,nx2) !ale
           CLOSE(1)
      ENDIF



!     End iteration
      END DO


      DEALLOCATE (wk,wh,k2,kx,ky,eku)

      STOP

      END PROGRAM hd2d
!=====================================================================
 
      SUBROUTINE evaluate(wk,dw,kx,ky,k2,dt,iadd)

      USE size_params
      USE hd_params
      USE fftw_plans
      IMPLICIT NONE
      INCLUDE 'fftw_include.h'
 
      COMPLEX, Intent(in)    ::  wk(nx1,ny)
      COMPLEX, Intent(inout) ::  dw(nx1,ny)
      REAL, Intent(in)       ::  kx(nx1), ky(ny), k2(nx1,ny)
      INTEGER, Intent(in)    ::  iadd
      REAL, Intent(in)       ::  dt
      
!     Local auxiliary arrays
      COMPLEX     ::  s1(nx1,ny),s2(nx1,ny),s3(nx1,ny)

      REAL        ::  dtnu,dtf0,dtnl
      INTEGER     ::  i,j

      dtnu=dt*xnu
      dtf0=dt*f0
      dtnl=dt*units

      if (iadd.eq.1) then
        dw=wk
      endif

!     Dissipative terms
      dw = dw - dtnu*k2*wk
     
!     Forcing
      do j=1,ny
       do i=1,nx1
         if ((k2(i,j).lt.kf02).and.(k2(i,j).gt.kf01)) then
            dw(i,j)=dw(i,j)+cmplx(dtf0,dtf0)
         endif
       end do
      end do

!     Non-linear terms
!     div (w {\bf v})
      CALL psik (wk,s2,k2)                                    ! psi(k)
      CALL nipgrd2D (s2,2,s1,kx,ky)                           ! vx(k)
      CALL rfftwnd_f77_complex_to_real(plan2DCR,1,s1,  &
                                  1,n2dh,0,0,0)               ! vx(x)
 
      s2 = wk                                                 ! w(k)
      CALL rfftwnd_f77_complex_to_real(plan2DCR,1,s2,  &
                                  1,n2dh,0,0,0)               ! w(x)

      call nipprod(s1,s2,s3)                                  ! vx(x) w(x)

      CALL rfftwnd_f77_real_to_complex(plan2DRC,1,s3,  &
                                  1,n2d,0,0,0)                !
      CALL ipgrd2D (s3,1,kx,ky)                               ! d/dx (vx w)

      dw = dw - dtnl * s3                                     ! dw update

                                                              ! s3 FREE
      CALL psik (wk,s3,k2)                                    ! psi(k)
      CALL ipgrd2D (s3,1,kx,ky)                               ! -vy(k)
      CALL rfftwnd_f77_complex_to_real(plan2DCR,1,s3,  &
                                  1,n2dh,0,0,0)               ! - vy(x)

      call ipprod (s2,s3)                                     ! -vy(x) w(x)

      CALL rfftwnd_f77_real_to_complex(plan2DRC,1,s2,  &
                                  1,n2d,0,0,0)                ! -(vy w)_k

      CALL ipgrd2D (s2,2,kx,ky)                               ! d/dy (- w vy)

      dw = dw + dtnl * s2                               ! NL w term COMPLETE !


      RETURN
      END SUBROUTINE evaluate
!=====================================================================
	SUBROUTINE globals(wk,k2,time)

	USE size_params
	USE hd_params
	IMPLICIT NONE

	COMPLEX, Intent(in)      ::  wk(nx1,ny)
	REAL, Intent(in)         ::  k2(nx1,ny)
	REAL, Intent(in)         ::  time

	REAL          ::  eneu,ens,pals,fct
	REAL          ::  wr,wi,ks
	INTEGER       ::  i,j,ny2

	!eneu = total kinetic energy = 1/2 sum (wk^2 / k^2)
	!ens  = total enstrophy = 1/2 sum (wk^2)
	!pals = total palinstrophy = 1/2 sum (k^2 * wk^2)
	!fct  = total forcing term = sum (f0k * wk / k^2)

	eneu   =  0.0
	ens    =  0.0
	pals   =  0.0
	fct    =  0.0

	ny2 = ny/2
	do j=1,ny
		do i=1,nx1
			!Count only the energy of the mode if ( kx > 0 .or. ky > 0 )
			if ( i.gt.1 .or. (j.gt.1 .and. j.le.ny2) ) then
				wr = real(wk(i,j))
				wi = aimag(wk(i,j))
				ks = k2(i,j)
				eneu = eneu + (wr*wr+wi*wi)/ks
				ens  = ens  + (wr*wr+wi*wi)
				pals = pals + (wr*wr+wi*wi)*ks
				if ((ks.gt.kf01).and.(ks.lt.kf02)) &
				fct  = fct  + (f0*wr+f0*wi)/ks
			endif
		end do
	end do
            
	!Write results in output file
	open(2,file='globals.out',position='append')
	!write(*,advance='no') '*'
	!write(*,'(A,$)') '*'
	write(2,*) time,eneu,ens,pals,fct
	close(2)

	RETURN
	END SUBROUTINE globals
!=====================================================================

      SUBROUTINE spec2D (ww,k2,eku)

      USE size_params
      IMPLICIT NONE

      COMPLEX, Intent(in)           ::   ww(nx1,ny)
      REAL, Intent(in)              ::   k2(nx1,ny)
      REAL, Intent(out)             ::   eku(nx2)
      REAL                          ::   wr,wi,ks,sk
      INTEGER                       ::   i,j,ny2,k

      eku = 0.

      ny2 = ny/2
      do j=1,ny
         do i=1,nx1
!        Count only the energy of the mode if ( kx > 0 .or. ky > 0 )
         if ( i.gt.1 .or. (j.gt.1 .and. j.le.ny2) ) then
            wr = real(ww(i,j))
            wi = aimag(ww(i,j))
            ks = k2(i,j)
            sk = sqrt(ks)
            if (sk .le. float(nx2)) then
               k = ifix(sk+.5)
               eku(k) = eku(k) + (wr*wr+wi*wi)/ks
            endif
         endif
         end do
      end do

      RETURN
      END SUBROUTINE spec2D
!=====================================================================

      SUBROUTINE init_kays (kx, ky, k2)

!-------------------------------------------------------
!     kx = 0,1,2,...,nx/2           
!     ky = 0,1,2,...,ny/2,-(ny/2-1),...,-2,-1    
!     k2 = kx*kx + ky*ky
!
!     Adjusted so that largest magn ky is POSITIVE
!-------------------------------------------------------

      USE size_params

      REAL, Intent( out), DIMENSION (nx1)        ::  kx
      REAL, Intent( out), DIMENSION (ny)         ::  ky
      REAL, Intent( out), DIMENSION (nx1,ny)     ::  k2

      INTEGER ::  i, j

      do i = 1, nx1
         kx (i) =  float(i-1)
      end do

      do j = 1, ny
         ky (j) =  float( MOD(j-1+ny/2,ny) - ny/2)
      end do
      ky (ny/2+1) =  float(ny/2)                       !.. Make ky_max > 0

      do j = 1, ny
        do i = 1, nx1
           k2 (i,j) = kx(i)*kx(i) + ky(j)*ky(j)
        end do
      end do

      END SUBROUTINE init_kays
!====================================================================

      SUBROUTINE psik (wk,fk,k2)

      USE size_params
      COMPLEX, Intent(in)      ::  wk(nx1,ny)
      COMPLEX, Intent(out)     ::  fk(nx1,ny)
      REAL, Intent(in)         ::  k2(nx1,ny)

      INTEGER                  ::  i, j

      do j=1,ny
         do i=1,nx1
            if (k2(i,j).ne.0) fk(i,j)=wk(i,j)/k2(i,j)
         end do
      end do
      fk(1,1)=0.0

      RETURN
      END SUBROUTINE psik
!====================================================================

      SUBROUTINE ipgrd2D (aa,idir,kx,ky)

      USE size_params

      COMPLEX, Intent(inout)   ::  aa(nx1,ny)
      REAL, Intent(in)         ::  kx(nx1),ky(ny)
      INTEGER                  ::  idir
      COMPLEX                  ::  ax

      if (idir.eq.1) then

      do j = 1, ny
         do i = 1, nx1
            ax = cmplx(0,-kx(i)) * aa(i,j)
            aa(i,j) = ax
         end do
      end do

      else 

      do j = 1, ny
         do i = 1, nx1
            ax = cmplx(0,-ky(j)) * aa(i,j)
            aa(i,j) = ax
         end do
      end do

      endif

      RETURN
      END SUBROUTINE ipgrd2D
!====================================================================

      SUBROUTINE nipgrd2D (aa,idir,ax,kx,ky)

      USE size_params

      COMPLEX, Intent(in)      ::  aa(nx1,ny)
      COMPLEX, Intent(out)     ::  ax(nx1,ny)
      REAL, Intent(in)         ::  kx(nx1),ky(ny)
      INTEGER                  ::  idir

      if (idir.eq.1) then

      do j = 1, ny
         do i = 1, nx1
            ax(i,j) = cmplx(0,-kx(i)) * aa(i,j)
         end do
      end do

      else 

      do j = 1, ny
         do i = 1, nx1
            ax(i,j) = cmplx(0,-ky(j)) * aa(i,j)
         end do
      end do

      endif

      RETURN
      END
!====================================================================
      
      SUBROUTINE ipprod (aprod, aa)

      USE size_params
      IMPLICIT NONE
      REAL, Intent(inout)        ::   aprod(2*nx1,ny)
      REAL, Intent(in)           ::   aa(2*nx1,ny)
      INTEGER                    ::   i,j

      do j=1,ny
         do i=1,2*nx1
            aprod(i,j) = aprod(i,j) * aa(i,j)
         end do
      end do

      RETURN
      END SUBROUTINE ipprod
!====================================================================
      
      SUBROUTINE nipprod (aa,ww,aprod)
      
      USE size_params
      IMPLICIT NONE
      REAL, Intent(in)           ::   aa(2*nx1,ny), ww(2*nx1,ny)
      REAL, Intent(out)          ::   aprod(2*nx1,ny)

!     Computes  aprod = aa(x) * ww(x), leaving aa and ww unchanged
!     Arrays are COMPLEX in main program, but REAL in this subr

      aprod = aa * ww

      RETURN
      END SUBROUTINE nipprod
!====================================================================

      SUBROUTINE preal2D (aa)
!     for use with (kx,ky) MHD2D code 
!
!     Imposes the interior AND edge reality conditions for a
!     pseudospectral algorithm.

      USE size_params
      IMPLICIT NONE

      COMPLEX,    Intent( inout)  :: aa(nx1,ny)
      INTEGER                     :: jy, jcjg
! ----------------------------------------------------

!     Modes (0, ky) are conj to (0, -ky)
      DO  jy = 2, ny/2
          jcjg = ny - jy + 2
          aa (1,jcjg) = CONJG ( aa (1,jy) )
      ENDDO

!     Do the special point k=(0,0) explicitly [=0, not just pure real]:
      aa(1,1) =  0.0

!     Zero faces
      aa (:, ny/2+1)  =  0.0
      aa (nx1, :) =  0.0


      RETURN
      END subroutine preal2D
!=====================================================================

      SUBROUTINE circ2D (ww,k2,k2lim)

      USE size_params
      IMPLICIT NONE

      COMPLEX, Intent(inout)     ::  ww(nx1,ny)
      REAL, Intent(in)           ::  k2(nx1,ny)
      REAL, Intent(in)           ::  k2lim
      INTEGER                    ::  i,j

      do j=1,ny
         do i=1,nx1
            if (k2(i,j) .gt. k2lim) then 
                ww(i,j)=0.0
            endif
         end do
      end do

      RETURN
      END SUBROUTINE circ2D
!=====================================================================
        
      SUBROUTINE set_plan1D (ny)

!     Set the FFTW 1D plans
!     Plans are passed to calling programs through the MODULE fftw_plans

      USE fftw_plans
      INCLUDE 'fftw_include.h'
      INTEGER, Intent( in)     ::  ny

      CALL fftw_f77_create_plan( plan1Df, ny, FFTW_forward,    &
                                FFTW_measure+FFTW_in_place)

      CALL fftw_f77_create_plan( plan1Db, ny, FFTW_backward,   &
                                FFTW_measure+FFTW_in_place)


      END SUBROUTINE set_plan1D
!=====================================================================

      SUBROUTINE set_plan2D (nx,ny)

!     Set the FFTW 2D plans
!     Plans are passed to calling programs through the MODULE fftw_plans

      USE fftw_plans
      INCLUDE 'fftw_include.h'
      INTEGER, Intent( in)  ::  nx, ny

      CALL rfftw2d_f77_create_plan( plan2DRC, nx, ny,       &
                 FFTW_real_to_complex, FFTW_measure+FFTW_in_place)

      CALL rfftw2d_f77_create_plan( plan2DCR, nx, ny,       &
                 FFTW_complex_to_real, FFTW_measure+FFTW_in_place)


      END SUBROUTINE set_plan2D
!=====================================================================

      SUBROUTINE write_xarr (aa,time)

      USE size_params
      IMPLICIT NONE
      REAL, Intent(in)           ::   aa(2*nx1,ny)
      REAL, Intent(in)           ::   time
      INTEGER                    ::   i,j

		!write(1) time,((aa(i,j),i=1,nx),j=1,ny) !ale
		DO i=1,nx
			DO j=1, ny
				write(1,'(F8.3,X$)') aa(i,j)
			end do
			write(1,*) 
		end do
		print*, time
		!ale

      RETURN
      END SUBROUTINE write_xarr
!=====================================================================

      SUBROUTINE read_xarr (aa,time)

      USE size_params
      IMPLICIT NONE

      REAL                        ::   aa(2*nx1,ny)
      REAL                        ::   time
      INTEGER                     ::   i,j

      read(1) time,((aa(i,j),i=1,nx),j=1,ny)
      aa(nx+1,:)=0.
      aa(nx+2,:)=0.

      RETURN
      END SUBROUTINE read_xarr
!=====================================================================

      SUBROUTINE initial_state_real (wk,u0)

      USE size_params
      IMPLICIT NONE

      REAL                        ::   wk(2*nx1,ny)
      REAL                        ::   u0
      REAL*8                      ::   dpi
      INTEGER                     ::   i,j

      dpi=4.d0*datan(1.d0)
      DO j=1, ny
         DO i=1,nx
            if (j .gt. ny/2) then 
               wk(i,j)=u0*dpi*((j-1)*2./float(ny)-1.)
            else 
               wk(i,j)=-u0*dpi*((j-1)*2./float(ny)-1.)
            endif
         END DO
      END DO

      wk(nx+1,:)=0.
      wk(nx+2,:)=0.

      RETURN
      END SUBROUTINE initial_state_real
!=====================================================================
